home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Pascal Super Library
/
Pascal Super Library (CW International)(1997).bin
/
MATH
/
PMAT12
/
PMAT.EXE
/
PMAT.DOC
< prev
next >
Wrap
Text File
|
1993-01-24
|
49KB
|
1,311 lines
P-Mat v1.2
An Turbo Pascal program for Recursive Matrix Algebra
Mark Von Tress, Ph.D.
PO Box 171173
Arlington TX 76003
Compuserve User ID : 71530,1170
Date: January 30, 1993
You are responsible for what you do with the
code. Here is a formal disclaimer:
DISCLAIMER: THIS PROGRAM IS PROVIDED AS IS, WITHOUT ANY
WARRANTY, EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED
TO FITNESS FOR A PARTICULAR PURPOSE. THE AUTHOR DISCLAIMS
ALL LIABILITY FOR DIRECT OR CONSEQUENTIAL DAMAGES RESULTING
FROM USE OF THIS PROGRAM.
Copyright (c) Mark Von Tress 1993
Mat-P 3
Table of Contents
I. INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . 4
II. FINANCIAL CONDITIONS OF USE . . . . . . . . . . . . . . . 4
III. BASIC ALGORITHMS . . . . . . . . . . . . . . . . . . . . 5
A. Structures and Allocation . . . . . . . . . . . . . . 5
B. The global stack: dispatch . . . . . . . . . . . . . 7
C. Recursion . . . . . . . . . . . . . . . . . . . . . . 7
D. Matrix Assignment . . . . . . . . . . . . . . . . . . 8
E. Parameter Passing . . . . . . . . . . . . . . . . . . 9
IV. FUNCTIONS . . . . . . . . . . . . . . . . . . . . . . . . 10
A. Error Handling . . . . . . . . . . . . . . . . . . . 10
B. Array Allocation or Deallocation . . . . . . . . . . 11
C. Stack Control . . . . . . . . . . . . . . . . . . . . 11
D. Matrix Allocation, Deallocation and Copy . . . . . . 13
E. Display Functions . . . . . . . . . . . . . . . . . . 14
F. Matrix IO . . . . . . . . . . . . . . . . . . . . . . 15
G. Binary Matrix Functions . . . . . . . . . . . . . . . 16
H. Unary Matrix Functions . . . . . . . . . . . . . . . 17
I. Patterned Matrices . . . . . . . . . . . . . . . . . 19
J. Mathematical Functions . . . . . . . . . . . . . . . 20
V. COMPILATION AND LIMITATIONS . . . . . . . . . . . . . . . 24
VI. REVISION HISTORY . . . . . . . . . . . . . . . . . . . . 25
A. Version 1.0 . . . . . . . . . . . . . . . . . . . . 25
B. Version 1.1 . . . . . . . . . . . . . . . . . . . . 25
C. Version 1.2 . . . . . . . . . . . . . . . . . . . . 25
Mat-P 4
I. INTRODUCTION
PMAT is Turbo Pascal (TP) source code for recursive matrix
algebra. The program is based on another program I wrote in C to
do the same thing. Both use a stack of matrices to keep track of
intermediate heap allocations. Once I figured it out in C it was
pretty easy to step back and do it in Pascal. The object
oriented extensions in TP also helped smooth the process. This
allowed a convenient scheme to allow matrices to be larger than
64K.
Of course you can't overload operators in TP, so the recursion is
messy. However, you can keep track of intermediate allocations on
the heap using PMAT.
new( a, makematrix( 1, 1 ) );
x = matequals(x, add( a,add(b,c)));
This code segment has to keep track of matrix allocations on the
heap, and then delete the temporary matrices. In this example,
the sum of B and C is a temporary matrix which would be lost
without some sort of global memory allocation tracking such as a
stack. Its memory should be deleted after the equals function is
called. The stack helps avoid the assignment of temporaries in
recursive calls.
On a more personal note, I wrote this program for fun. I started
on this problem in TP 3.0 when that was the best compiler on the
market (1984). I never really got a satisfactory answer, so I set
it down. Then I picked it up in Turbo C, and didn't get a
satisfactory answer. Then I picked it up in C++, and got a good
answer in about 1990 (see YAMP on the BPROGB of Compuserve). Then
I wrote it in C with good results, then TP 6.0 just for old times
sake. I had forgotten how fast TP compiles. It's a real code
blaster. Object oriented programming in TP 6 also helped. I guess
I learned some computer science along the way too.
II. FINANCIAL CONDITIONS OF USE
If you like this program, send me $5.00 US (or buy a deserving
beggar lunch. In either case, find a way to repay my charity.)
You may distribute the package unchanged. I only plan to
distribute it electronically. I retain the copyright as proof of
authorship.
Mat-P 5
III. BASIC ALGORITHMS
A. Structures and Allocation
PMAT has two important objects: vmatrix and mstack. The vmatrix
object is declared in the interface
vmatrixptr = ^vmatrix;
vmatrix = Object
r, c : integer;
Function m( i, j: integer ): double;
Function mm( i, j: integer ) : fp;
constructor MakeMatrix( vr, vc: integer );
Destructor KillVmatrix;
Procedure Garbage;
Procedure Show( strng: String );
Procedure infomatrix( strng: String );
private
v,vcheck: ^app;
nelems : longint;
onstack : boolean;
Procedure purgevectors;
Procedure allocvectors( rr, cc: integer );
End;
The vmatrix contains integers for the dimensions of the matrices,
an array of pointers to vectors of doubles, the number of
elements, and a check address. Matrix allocation is based on
Numerical Recipes in C . An array of addresses of vectors of
doubles is allocated first. Then, an array is allocated and its
first address is stored in the address vector. This method allows
matrices larger than 64K and scatters the rows of the matrix
throughout the heap. A row is restricted to 8190 doubles to keep
it under 64K.
Object oriented programming helped hide the ugly notation
required for storing and retrieving elements from a vmatrix. You
use the member function mm(i,j) to store elements in the matrix,
and m(i,j) to retrieve elements. An example is
x^.mm(i,j)^ := y^.m(i,j);
where x and y are vmatrix pointers. Note the extra ^ in
x^.mm(i,j)^. This returns the address in the heap for storing an
element. Using x^.mm(i,j) won't work. Think of element assignment
as saying "put y^.m(i,j) in the address pointed to by
x^.mm(i,j)". Using m(i,j) returns a double via
Mat-P 6
m := v^[i]^[j];
and using mm(i,j) returns an address via
mm := @v^[i]^[j];
Using member functions protects you from having to type these
unintuitive indexing operations. I also allows easy access to
matrices larger than 64K. They also do range checking on the
indexes.
The check pointer contains the value of v upon allocating the
matrix. Then it is set to zero upon deallocation of v. The check
value is used in the member function Garbage(). If v does not
equal vcheck, then the matrix has been deallocated, or memory may
have been corrupted. The program halts if you try to use a
garbage matrix (better safe than sorry!).
The boolean variable, onstack, indicates that the matrix is on
the stack (onstack = true) or not (onstack = false). This
controls how matrices are copied. Its initial value is false. It
is set to true by push(), and false by pop(). See CopySD().
PMAT expects you to use pointers to matrices, so allocation looks
like this
procedure junk;
var x,y,z : vmatrixptr;
begin
new( x, makematrix( 1, 1 ) );
new( y, makematrix( 1, 1 ) );
new( z, makematrix( 1, 1 ) );
.
.
.
{ put code that uses x,y,z }
.
.
.
dispose( x, killvmatrix );
dispose( y, killvmatrix );
dispose( z, killvmatrix );
end;
Call the constructor for each matrix, then make sure to dispose
of the storage at the end of the procedure. Your program will
develop memory leaks if you do not dispose of the vmatrix
pointers. PMAT is not consistent with the MARK-RELEASE method of
allocating memory.
Mat-P 7
The mstack type is
mstackptr = ^mstack;
mstack = Object
constructor InitMatrixStack;
Destructor KillMatrixStack;
Procedure inclevel;
Procedure declevel;
Function ReturnMat: vmatrixptr;
Function DecReturn: vmatrixptr;
Procedure push( Var a: vmatrixptr );
Procedure nrerror( strng: String );
Procedure DumpStack;
private
nummats, level : integer;
next : mstackptr;
mv : vmatrixptr;
Procedure cleanstack( a: vmatrixptr );
Function pop: vmatrixptr;
End;
The mstack structure is a typical single linked list sort of
thing. The nummats field contains the number of matrices on the
stack. The level field contains the subroutine nesting level
which keeps track of which temporary matrices to be disposed. As
you enter and leave functions that use matrix assignment, you
increment and decrement the function nesting level. The
mstackptr, next, is a pointer to the next mstack object. The
vmatrixptr mv is the address of a vmatrix object referenced by
the stack object.
B. The global stack: dispatch
Mstack objects are stored in the global stack linked list with a
base pointer called dispatch. Dispatch is automatically
initialized in the use statement.
Dereference items in dispatch using the ^. operator, such as
dispatch^.level. The stack also has a tail, so if you try to
deallocate it, the program stops. It has a value of nil in the
field called next.
C. Recursion
The recursive matrix functions push mstack pointers onto the
stack. The address of the matrix reference is stored in mv. You
create stack elements when you call the push function. mstack
elements are deallocated automatically by other functions, so you
Mat-P 8
don't have to do it. Popping occurs during calls to matequals()
for matrix assignment. Temporary mstack pointers are freed during
matequals by the garbage collector function, CleanStack(). You do
not have to pop explicitly.
The only part of mstack that you have to pay attention to is the
function nesting level. Garbage collection occurs during
assignment. Any function that uses matrix assignment pushes
matrices onto dispatch. Stack elements are deleted during garbage
collection if their value of level is greater than or equal to
the value of level held in dispatch. This way, each function has
its own stack, but all stack elements are stored in the same
linked list.
You control level with the member functions, Inclevel(),
Declevel(), and DecReturn(). Inclevel() and Declevel() increment
or decrement level. Declevel() checks that level has not become
negative. The program stops if it has.
DecReturn() is used for returning matrices from functions. It
decrements the nesting level, and returns the address of the
matrix, v, in dispatch^.next. An example of returning a matrix
pointer recursively is
function f1: vmatrixptr;
var x: vmatrixptr;
begin
dispatch^.inclevel;
new( x, makematrix(1,1));
x := matequals(x, ident(5) );
dispatch^.push(x);
f1 := dispatch^.decreturn;
end;
The stack works by manipulating pointers to vmatrix objects,
rather than copies. This is faster and uses less memory. The
matrix pointer is pushed onto the stack. Then the function level
is decremented, and the pointer x is returned as
dispatch^.next^.mv. There is no need to dispose x, since it is
just a pointer stored in a global stack. The stack functions will
dispose it when it is appropriate.
D. Matrix Assignment
Matrix assignment is accomplished using the function matequals,
having prototype
Function matequals(Var lop: vmatrixptr; rop: vmatrixptr) :
vmatrixptr;
The statement from the example above uses assignment
Mat-P 9
x := matequals( x, ident(5) );
The matrix being assigned must be used in two places. The first
place is in the function call where the contents of the matrix
objects are manipulated. The second is on the left of the equals
sign where the current function is told where the new address of
x is to be found.
Matequals first checks that the two inputs are not garbage. Then
it checks if the two inputs point to the same 2 dimensional
array. If they do, then the matrix a gets a new 2 dimensional
array. The contents of rop are copied into lop. Then the
intermediate stack matrices are deleted.
E. Parameter Passing
Parameter passing of matrices is the same as in ordinary Pascal.
If you want to change the contents of a matrix that was passed as
a parameter, then declare it a var parameter in the procedure.
Make sure that the matrix is correctly initialized before the
call. Then use matrix assignment to make the changes in the
matrix. Here is an example.
procedure f1( var a: vmatrixptr );
begin
dispatch^.inclevel;
a := matequals(a,ident(3));
dispatch^.declevel;
end;
procedure f2( void );
var x:vmatrixptr;
begin
dispatch^.inclevel;
new(x,makematrix(1,1));
x^.show('before pass');
f1(x);
x^.show('after pass');
dispose(x,killvmatrix);
dispatch^.declevel;
end;
The call to f1 in f2() passes x as a var vmatrix pointer. The
function call changes the contents of x from a 1x1 matrix to the
3x3 identity matrix.
This example clarifies why matequals() needs to use the formal
parameter name on the left of the equals sign. You have to tell
the program where you put the new copy of the vmatrix.
Mat-P 10
IV. FUNCTIONS
There are several kinds of functions in PMAT. They fall in
several categories:
0. error handling
1. array allocation or deallocation
2. stack control
3. matrix allocation, deallocation and copy
4. display
5. IO
6. binary matrix functions
7. unary matrix function
8. patterned matrices
9. mathematical functions
The interface should provide a quick reference to the names and
spellings. The following discussion explains the header.
A. Error Handling
Procedure nrerror( strng: String );
This is the error handler from Numerical Recipes in C. It is
called whenever a fatal error occurs. It prints the error text
and exits with a value 1. It is a member function of mstack
objects. Call it as dispatch^.nrerror('something''s happening
here');
Procedure Garbage;
This function checks if the matrix a has been deleted or
corrupted. It calls nrerror if it has. This is a vmatrix member
function. Call it as x^.Garbage;
Mat-P 11
B. Array Allocation or Deallocation
The declarations form the structure for matrices.
ap = Array[1..1] Of double;
apptr = ^ap;
app = Array[1..1] Of apptr;
The type ap is an array of doubles. Note that the range is from 1
to 1, and the array is subscriptable and variable length. Range
checking is turned off in the matrix allocation and indexing
functions so that these arrays act like they do in C. Each row of
the matrix is allocated separately by new(apptr). Each new apptr
is stored in an subscripatble, variable length array of pointers
to arrays of doubles.
Procedure allocvectors( rr, cc: integer )
Arrays are allocated using the Numerical Recipes in C method for
dmatrix arrays. Matrices in PMat are indexed from 1 to r, and 1
to c.
procedure purgevectors.
Arrays are freed using the Numerical Recipes in C method for
dmatrix arrays.
Neither of these functions need to be called directly since they
are called by the vmatrix construction routines. They are private
to the vmatrix class, so they cannot be used outside of the
pmat.pas unit.
C. Stack Control
dispatch: mstackptr; (* stack top: global *)
This declaration in the interface makes dispatch available to all
units that use pmat. It is allocated automatically from the use
statement.
constructor InitMatrixStack; (* set up stack top and tail *)
This function allocates dispatch and the stack tail. It is called
by the pmat initialization routine.
Procedure push( Var a: vmatrixptr ); { push a matrix onto stack }
pushes a matrix onto the stack. It should be called before a f1
Mat-P 12
:= dispatch^.ReturnMat; or f1 := dispatch^.DecReturn; call. The
new matrix is stored in dispatch^.next^.mv. push() also sets
a^.onstack to true. This is a member function of mstack.
Function pop: vmatrixptr; (* free dispatch->next, return v*)
pop is a private function of the matrix stack. You should not
call it. It is called repetitively by CleanStack to delete
temporary matrices. It deletes the matrix structure, and returns
the matrix pointer mv. pop() also sets mv^.onstack to false.
vmatrix *ReturnMat( void ); /* return stack top matrix */
ReturnMat() returns the value dispatch^.next^.v in functions that
have not incremented the function nesting level, i.e. required
matrix assignment. It does not modify the subroutine nesting
level. All of the binary matrix functions push a vmatrix pointer
onto the stack, and then return by the statement "f1 :=
dispatch^.ReturnMat;". It is a public member function of mstack.
Function DecReturn: vmatrixptr; (* return stack to
matrix,decrement level*)
DecReturn returns the value dispatch^.next^.v in functions that
have incremented the function nesting level, i.e. required matrix
assignment. It decrements the subroutine nesting level. Functions
that use matrix assignment, such as inv(), push matrix pointers
onto the stack, and then return by the statement
"f1 := dispatch^.DecReturn;". This is a public member function of
mstack.
Procedure inclevel; (* increment function nesting level *)
Inclevel() increments the function nesting level. See the
discussion of the stack and recursion. This a public member
function of mstack, and is called as dispatch^.inclevel.
Procedure declevel; (* decrement function nesting level *)
Declevel() decrements the function nesting level. It also checks
that the level has become negative. If so, it calls nrerror().
This indicates that level has been decremented more often than it
has been incremented. Check the Inclevel()-Declevel() pairs, or
the Inclevel()-DecReturn() pairs. See the discussion of the stack
and recursion. This a public member function of mstack, and is
called as dispatch^.declevel. Use this function when you use
Mat-P 13
matrix assignment in a procedure, but do not want to return a
matrix.
Procedure cleanstack( a: vmatrixptr );
(* pop stack elements if level >= dispatch^.level *)
CleanStack() is a private function of the stack. It should not be
called by the user. The argument is a matrix that is not to be
freed. Thus the call "CleanStack( x )" indicates that all
matrices with level at or above dispatch^.level should be freed,
except for matrices that have the same address as x. It calls pop
while (m^.level >= dispatch^.level), and (dispatch^.next^.next !=
NIL). It frees all matrices not equal to x.
Function matequals( Var lop: vmatrixptr; rop: vmatrixptr ) :
vmatrixptr; (* copy rop into lop *)
The matrix being assigned must be used in two places. The first
place is in the function call where the contents of the matrix
structure are manipulated. The second is on the left of the
equals sign where the current function is told where the new
address of x is to be found.
matequals() first checks that the two inputs are not garbage. The
contents of rop are copied into lop by CopySD(). Then the
intermediate stack matrices are deleted by CleanStack().
D. Matrix Allocation, Deallocation and Copy
constructor MakeMatrix( vr, vc: integer );
(* allocate a matrix *)
MakeMatrix() constructs a matrix with vr rows and vc columns.
This should be called to properly construct a matrix before using
it in the matrix functions. This is a public member function of
vmatrix, and should be called as new(x, makematrix(1,1));
Destructor KillVmatrix; (* deallocate a matrix *)
This frees the matrix storage in m^.v, and then dispose m. It
calls nrerror if m is corrupted. It is a public member function
of vmatrix and should be called as dispose(m,killvmatrix);
Procedure CopySD( Var source, dest: vmatrixptr );
(* copy a matrix from source to dest *)
This function should be considered private to the stack function,
Mat-P 14
matequals(). It is also hidden to the user in the implementation
of pmat.
CopySD() replaces the contents in the destination matrix by a
copy of the contents in the source matrix. If source and dest are
the same, then CopySD() returns immediately since you are trying
to copy a matrix onto itself.
CopySD() recognizes when matrices are being copied from the
stack. Matrices being copied from the stack are about to be
deleted by Cleanstack(), so their matrix pointer is popped into
the destination matrix pointer. The destination matrix storage is
released first. The indices are reset. Next, the stack matrix
pointer is popped into the destination matrix pointer. Then the
function returns.
If the matrices point to the same dmatrix, and the source matrix
is not on the stack, then a new dmatrix array is allocated for
dest. If they have a different number of rows or columns, the
destination 2D array is reallocated to have the same number of
rows and columns. Then the dmatrix array contents of the source
are copied into the destination.
E. Display Functions
Procedure DumpStack;
This dumps the stack. It is the only function that may be called
between push and the function return. It is for debugging
purposes. It is a public member function of mstack and should be
called as dispatch^.dumpstack;
Procedure infomatrix( strng: String ); (* show dimensions *)
infomatrix() displays information about the matrix. The strng is
a string describing the matrix. The string, strng, is displayed,
then the number of rows, columns, m, and mcheck. Use
infomatrix() when the matrix is large, or when you are having
problems with corrupted or deleted matrices. This is a public
function of vmatrix and should be called as x^.Infomat('something
is happening here');
Mat-P 15
Procedure Show( strng: String );(* display matrix *)
show() displays a matrix. The string, strng is also displayed
before the matrix elements. You can control the width and decimal
places by calling setwid() and setdec() before calling show().
This is public function of vmatrix and a typical call is
"x^.show("x");"
Procedure setwid( wid: integer ); (* set display width *)
Procedure setdec( decimals: integer ); (* set display decimals *)
Function getwid: integer; (* get display width *)
Function getdec: integer; (* get display decimals *)
setwid() and setdec() store hidden global integers for the
display width and decimals. show() and Writea() call getwid() and
getdec() to use these variables.
F. Matrix IO
External matrices are stored in an ASCII character format. The
first row is a title string that is at most 80 characters
including the final carriage return, '\n'. The second row
contains two integers separated by white space. The first integer
is the number of rows, r, and the second is the number of
columns, c. The next rows are the rows of the matrix. Matrix may
span several text rows, but there must be r*c elements.
Function reada( fid: String ): vmatrixptr;
(* read ascii matrix *)
Reada() reads a matrix stored in the ASCII file pointed to by
fid. Reada() terminates if the title string is longer than 80
total characters, including the terminal carriage return, '\n'.
Reada() should be called in a matrix assignment since it returns
a matrix off of the stack. As an example,
x := matequals( x, Reada('xmatrix.dat');
Procedure writea( fid: String; a: vmatrixptr; titlestr: String );
(* write ascii matrix *)
Writea() stores a vmatrix, a, in file, fid, using the PMat ASCII
file format. Data is written to fid using the current print
display values found in setdec() and setwid(). The title string
Mat-P 16
is truncated to 80 characters (including '\n') if it is longer
than 80 characters.
G. Binary Matrix Functions
The binary matrix functions are recursive in that they push
matrices onto the stack, and do not manipulate the function
nesting level. They all push a vmatrix pointer onto the stack,
and return a vmatrix pointer using a call to ReturnMat(). The
incoming matrices are checked by calling Garbage() before doing
any calculations. A typical use of these functions is
(* d= a+(b-c) ; *)
d = equals( d, add(a,sub(b,c)));
The call to equals deletes the temporary difference of c from b.
Scalars may also be used in binary operations. Each of the binary
matrix functions have "sc" appended to the left(right) if a
scalar is to be added to a matrix on the left(right).
function add( a, b: vmatrixptr ): vmatrixptr; (* addition a+b *)
function scadd( a: double; b: vmatrixptr): vmatrixptr;
function addsc( b: vmatrixptr; a: double); vmatrixptr;
Matrix addition, conditions : a^.r == b^.r and a^.c == b^.c
The scalar additions add a double to a matrix.
Function sub( A, B: vmatrixptr ):vmatrixptr; (* subtraction a-b*)
Function scsub( a: double; b: vmatrixptr): vmatrixptr;
Function subsc( b: vmatrixptr; a: double); vmatrixptr;
Matrix subtraction, conditions : a^.r = b^.r and a^.c = b^.c
The scalar subtractions perform the subraction in the indicated
order: scsub(i,j) = A - b(i,j), subsc(i,j)= b(i,j)-A.
Function Mult( A, B: vmatrixptr ):vmatrixptr;(*matrix mult a*b*)
Caley multiplication, inner products are formed using extended
doubles, then truncated back to doubles when stored in a matrix.
conditions : a^.c = b^.r.
Mat-P 17
Function emult( a, b: vmatrixptr ):vmatrixptr;
(* elementwise mult a#b *)
function scmult( a: double; b: vmatrixptr): vmatrixptr;
function multsc( b: vmatrixptr; a: double); vmatrixptr;
Elementwise multiplication,
conditions : a^.r = b^.r and a^.c = b^.c
The scalar multiplications multiply the elements of b by a.
Function divis(a,b:vmatrixptr): vmatrixptr;
function scdivis( a: double; b: vmatrixptr): vmatrixptr;
function divissc( b: vmatrixptr; a: double); vmatrixptr;
Elementwise division,
conditions : a^.r = b^.r and a^.c = b^.c and for all i,j
b^.m(i,j) <> 0
The scalar divisions perform division in the indicated order:
scdivis(i,j) = A/b(i,j), divissc(i,j)= b(i,j)/A.
H. Unary Matrix Functions
The unary matrix functions, neg(), tran(), and inv(), are
recursive in that they push matrices onto the stack, and do not
manipulate the function nesting level. neg(), and tran() push a
vmatrix pointer onto the stack, and return a vmatrix pointer
using a call to ReturnMat(). inv() uses matrix assignment so it
increments the nesting level, pushes the return matrix, and
returns by a call to DecReturn(). The incoming matrices are
checked by calling Garbage() before doing any calculations. A
typical use of these functions is
Mat-P 18
function regression: vmatrixptr;
var x,y,xpx,xpy,beta : vmatrixptr;
begin
(* beta = inv(x'x)x'y *)
dispatch^.Inclevel;
new(x, makematrix(1,1));
new(y, makematrix(1,1));
new(xpx, makematrix(1,1));
new(xpy, makematrix(1,1));
new(beta, makematrix(1,1));
x := matequals( x, Reada( "x.dat") );
y := matequals( y, Reada( "y.dat") );
xpx := matequals( xpx, mult( tran(x), x ));
xpy := matequals( xpy, mult( tran(x), y ));
beta := matequals( beta, mult( inv(xpx), xpy));
dispose( x, killvmatrix);
dispose( y, killvmatrix);
dispose( xpx, killvmatrix);
dispose( xpy, killvmatrix);
dispatch^.push( beta );
regression := dispatch^.DecReturn;
end;
function sweepreg: vmatrixptr;
var xy, xypxy, beta : vmatrixptr;
begin
(* beta using sweep *)
dispatch^.Inclevel;
new(xy, makematrix(1,1));
new(xypxy, makematrix(1,1));
new(beta, makematrix(1,1));
xy := matequals( xy, Reada( "xy.dat") );
xypxy := matequals( xypxy, mult( tran(xy), xy ));
xypxy := equals( xypxy, sweep( xypxy, 1, xy^.c - 1) );
beta := matequals( beta,
submat( xypxy, 1, xy^.c -1, xy^.c, xy^.c) );
dispose( xy, killvmatrix);
dispose( xypxy, killvmatrix);
dispatch^.push( beta );
sweepreg := dispatch^.DecReturn;
end;
Function neg( A: vmatrixptr ): vmatrixptr;(* negation *)
Changes the sign of all elements.
Function tran( a: vmatrixptr ): vmatrixptr;(* transpose *)
Matrix transpose: c^.m(i,j) = a^.m(j,i)
Mat-P 19
Function sweep( A: vmatrixptr; k1, k2: integer ) : vmatrixptr;
Sweep cols k1...k2 along the diagonal and in place in a. See the
example for how to use sweep for regression. Input matrix must be
square. Row and columns having pivot elements less than 10E-8 are
set to zero. This indicates a colinearity with at least one
column that has already been swept out.
Function inv( Amat: vmatrixptr ): vmatrixptr;
(* inversion by sweep *)
Calls sweep for inversion: x = matequals( x, sweep(x, 1, x^.c ));
This is a Gaussian elimination inverter without column pivoting.
It fails to accurately invert Hilbert matrices of dimension 8 or
more, which is reasonable for inversion by Gaussian elimination.
I. Patterned Matrices
The patterned matrix functions are recursive. They push matrices
onto the stack and return by calling ReturnMat. They should be
used in conjunction with matrix assignment.
Function submat( a: vmatrixptr; lr, r, lc, c: integer ):
vmatrixptr; (*submatrix*)
Extract a submatrix of a, using rows lr to r, and columns lc to
c. All integer arguments must be in the range of the dimension of
a.
Function ident( n: integer ): vmatrixptr; (*identity matrix*)
Construct an n dimensional Identity matrix
Function fill( rr, cc: integer; d: double ):vmatrixptr;
(* matrix of constants d *)
Construct a rr x cc matrix of constants, d.
Function cv( a, b: vmatrixptr ): vmatrixptr;
(* vertical concatenation *)
Vertically concatenate b to a. a is on top, b is on bottom, and
they must have the same number of columns.
Mat-P 20
Function ch( a, b: vmatrixptr ): vmatrixptr;
(* horizontal concatenation *)
Horizontally concatenate b to a. a is on the left, b is on the
right, and they must have the same number of rows.
vmatrix *vecdiag( vmatrix *a );
(* place vector on diagonal *)
Place elements of a rx1 or 1xc matrix on the diagonal of a square
matrix.
J. Mathematical Functions
The mathematical function module for YAMP was translated for
PMAT. These are some of the many linear algebra functions that
are used.
function MSort(a :vmatrixptr; col, order: integer): vmatrixptr;
Sort the rows of a matrix by a column if order is not equal to
zero. If order is zero, then sort the columns of a matrix by a
row.
function Mexp(ROp: vmatrixptr): vmatrixptr;
function Mabs(ROp: vmatrixptr): vmatrixptr;
function Mlog(ROp :vmatrixptr): vmatrixptr;
function Msqrt(ROp: vmatrixptr): vmatrixptr;
These functions operate on the elements of a matrix. They halt if
an argument is out of range. The take exp(a^.m(i,j)),
abs(a^.m(i,j)), ln( a^.m(i,j)), and sqrt(a^.m(i,j)). You might
consider adding the trig functions too.
function Trace(ROp: vmatrixptr): double;
Takes the sum of the diagonal elements of a square matrix.
function Helm (n: integer): vmatrixptr;
Returns the Helmert matrix of dimension n.
function Index( lower, upper, step : integer): vmatrixptr;
Mat-P 21
Returns an index matrix with elements between lower and upper,
and in increments of step. If lower is less than upper, then step
must be positive. If lower is greater than upper, then step must
be negative.
function Kron(a,b:vmatrixptr): vmatrixptr;
Returns the Kroniker product of a and b. The blocks of the
Kroniker product are a^.m(i,j)*b.
Mat-P 22
function Det(Datain : vmatrixptr): double;
Returns the determinant of Datain, which must be a square matrix.
function Cholesky(ROp: vmatrixptr): vmatrixptr;
Returns the Cholesky decomposition of a square symmetric matrix:
ROp = S'S, where S is an upper triangular matrix. Cholesky stops
if ROp is singular. This routine does not use pivoting.
procedure QR(var ROp, Q, R: vmatrixptr; makeq: boolean);
This procedure produces the QR decomposition of an rxc matrix:
ROp=QR, where Q is an rxc matrix of orthogonal columns, and R is
an upper triangular cxc matrix. This routine stops if a zero is
detected along the diagonal of R. The variable makeq determines
if Q is produced. The procedure calculates Q = ROp*Inv(R). If the
calculation of Q is not required, or inversion of R is
prohibative, then set makeq = FALSE to avoid the inversion of R.
function Svd(var A, U, S, V : vmatrixptr;
makeu, makev: boolean) : integer;
This function returns an integer indicating that the singular
value decomposition has failed. The singular value decomposition
is calculated as A = UDiag(S)V' where U and V are orthogonal
matrices and S is a column vector containing the singular values
of A. The singular values of A are the eigenvalues of A'A. If you
wish to calculate U or V, then set makeu = TRUE or makev = TRUE
respectively.
function Ginv( ROp : vmatrixptr): vmatrixptr;
This function returns the generalized inverse from the singular
value decomposition of ROp=UDiag(S)V'. The generalized inverse of
ROp is Ginv(ROp) = VDiag(S+)U', where S+ is a diagonal matrix of
the reciprocals of S. If S has a element of zero, the
corresponding element of S+ is also zero.
function Fft(ROp : vmatrixptr; INZEE : boolean): vmatrixptr;
This function returns the fast Fourier transform of an input
matrix ROp. The length of ROp may be any value. ROp may have 1 or
2 columns. If it has 1 column, the vector is assumed to be a real
vector. If the matrix has 2 columns, the matrix is assumed to
consist of complex numbers. The first column contains the real
part, and the second column contains the imaginary part. The
Mat-P 23
forward transform is calculated if INZEE is TRUE, and the inverse
transform is calculated if INZEE is FALSE.
function Vec( ROp : vmatrixptr): vmatrixptr;
This functions stacks the columns of a matrix into a single
column vector. Thus if ROp is rxc, then Vec(ROp) has r*c rows and
1 column. The c blocks of Vec(ROp) are columns of ROp.
function Diag(ROp : vmatrixptr): vmatrixptr;
This function has two purposes. If the input matrix is a column
or row vector, then Diag places the elements along the diagonal
elements of square matrix of zeros. If the input matrix is
square, then all of the off-diagonal elements are set to zero.
function Shape(ROp: vmatrixptr; nrows: integer): vmatrixptr;
This functions reshapes a matrix to have nrows rows. If ROp is
has r*c elements, then nrows must divide r*c without a remainder.
type margins = (ALL,ROWS,COLUMNS);
This enumerated type controls the actions of the next 5
functions. If you use ALL, then the functions operate over all
elements of the input matrix. If you use ROWS, then the functions
operate over the rows of the input. If you use COLUMNS, then the
functions operate over the columns of the matrix.
function Sum( ROp : vmatrixptr; marg : margins ) : vmatrixptr;
This functions returns sum of the elements in a matrix. If
marg=ALL, then the returned matrix is 1x1, containing the sum of
all elements. If marg=ROWS then the return matrix has the same
number of columns as ROp, and 1 row. Each column contains the sum
of the row elements. If marg=COLUMNS then the return matrix has
the same number of rows as ROp, and 1 column. Each row cantains
the sum of the column elements.
function Sumsq( ROp : vmatrixptr; marg: margins): vmatrixptr;
This function is similar to Sum, except that it returns the sum
of squared elements.
Mat-P 24
function Cusum( ROp : vmatrixptr): vmatrixptr;
This function returns the cumulative sum across rows of matrix.
function Mmax( ROp: vmatrixptr; marg : margins): vmatrixptr;
function Mmin( ROp: vmatrixptr; marg : margins): vmatrixptr;
These functions find the max and min elements of a matrix, and
the position of the elements. If marg=ALL, then the functions
return a 3x1 matrix. The first element is the row index of the
max, the second element is the column index of the element, and
the third element is the max(or min). If marg=ROWS, then the
returned matrix is a 3xc matrix where the first row is the row
index of the max(min) element in a column, the second row is the
column index of the max(min) element, and the third row are the
max(min) for the column. Similarly, if marg=COLUMNS, then a rx3
matrix is returned.
Procedure Crow(var ROp: vmatrixptr; row : integer; c : double);
Procedure Srow(var ROp : vmatrixptr; row1, row2 : integer);
Procedure Lrow(var ROp : vmatrixptr; row1, row2 : integer; c :
double);
Procedure Ccol(var ROp : vmatrixptr; col : integer; c : double);
Procedure Scol(var ROp : vmatrixptr; col1, col2 : integer);
Procedure Lcol(var ROp : vmatrixptr; col1, col2 : integer; c :
double);
These functions perform elementary row and column operations on
matrices. They change the input matrices. The first letter
indicates the operation. The next 3 letters indicate if the
operation is for rows or columns. For Crow, the row indexed by
row is multiplied by c. For Srow, row1 and row2 are swapped. For
Lrow, row1 is multiplied by c, and added to row2. Ccol, Scol, and
Lcol perform similar operations on columns.
V. COMPILATION AND LIMITATIONS
Programs using PMAT must include pmat in the uses statement. Also
put pmatop in the uses statement if you use code from PMATOP.PAS.
Compile and link the program. See mattest2.pas or testop.pas for
example files using PMat or PMatop. Also see makefile for a
Borland make file for using the command line compiler.
Mat-P 25
PMat does have some limitations. It is limited by the available
memory in the machine, so large matrices will exhaust the heap,
and the program will stop. This is bothersome on a PC without a
DOS extender. The BP7 protected mode for DOS programs allows up
to 16 megs of memory, and matrices may be larger than 640K. Note
that a 1000x1000 of double takes 8 megs, so PMat may not be
suitable for really large matrices.
Another limitation is its readability. Using function calls
rather than overloaded algebraic operators muddles things. There
is no other obvious way to do this in TP though. One trick I use
to help me with this is to count open-close parenthesis pairs.
Some programming editors help this too by blink the other member
of the pair. Another trick is to keep recursive expressions on a
single line. This helps with storage also, since the stack
storage is cleaned up after each call to equals. If you have
memory overflow problems, try splitting the recursive expressions
into several parts by paying attention to where the big matrices
are created.
VI. REVISION HISTORY
A. Version 1.0
This is the original version. It is a copy of my matrix program
C-Mat. I may extend PMAT to have most of the capability of YAMP,
my C++ program. This would involve a virtual memory manager, a
graphics module, and larger set of mathematical functions. I
don't know if it is worth the effort though.
B. Version 1.1
This version adds BP 7.0 compatability for protected mode DOS
applications. The version 7 DPMI extender eliminates the need of
a virtual memory manager. You need a 286 or 386 to use this
feature.
C. Version 1.2
Added a new module for the math functions from YAMP. This
includes matrix sorts, elementwise functions, trace, helmert
matrix, index matrix, Kroniker products, determinants, Cholesky
decomposition, QR decomposition, singular valued decomposition,
generalized inverse, fast Fourier transform, vec, diagonal,
shape, sums, sums of squares, cumulative sums, max, min, and
elementary row-column operations. Added arithmetic operations
with scalars: scadd, addsc, etc.